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Abstract 

The parameter space of dynamical systems arising in applications is often found to be high- 
OO ' dimensional and difficult to explore. We construct a fast algorithm to numerically analyze "quan- 

fy^ ' titative features" of dynamical systems depending on parameters. Using a classical problem from 

^^ , mathematical ecology as an example, we demonstrate how to apply the algorithm to investigate 

t^ ■ the amplitude of a limit cycle depending on seven parameters. We stress the practical value of 

the algorithm but we also provide a rigorous error analysis to justify the overall strategy. Our 
approach turns out to be particularly useful in the case of comparing experimental data to a 
model defined by differential equations and to investigate whether the equations can approximate 
the modeled system. 

1 Introduction 
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^H ' In this paper we investigate the problem of how to explore parameter spaces for dynamical systems. 

r^ . We shall focus on ordinary differential equations (ODE), but remark that our approach works also in 

"ti ' different settings such as difference equations, partial differential equations or stochastic dynamics 



> 

m 



s 



Parametrized dynamical systems in mathematical modelling not only lead to questions about 
structural features of the model and their dependence on parameters, which is the is central question 
in applied bifurcation theory, but also to quantitative questions. The main question we are going 
to address is how a feature, like a limit cycle or a homoclinic orbit, depends on parameters. Some 
questions often encountered are: 

• What is the amplitude of a limit cycle? 



C^ ' • What is the time required to traverse a part of a homoclinic orbit? 

OO , • How far does an invariant manifold extend in phase space? 

All answers are obviously parameter dependent and are going to have relevant interpretations in the 
modelling problem under consideration. Drawing quantitative conclusions requires understanding 
^\ ' how different parameters influence the feature. Even if the parameter space is of moderate dimension 

(« 15 — 25 parameters) it is often infeasible to numerically explore it systematically by computing a 
dense grid of parameters values. In many mathematical models no dedicated techniques are used to 
explore the parameter space. Instead several "ad-hoc" methods are used which can potentially yield 
incorrect conclusions. 

Most existing techniques for exploring parameter spaces systematically are based on purely statis- 
tical methods and in particular on sensitivity analysis. A short exposition of the methods of sensitivity 
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Figure 1: Limit cycle at default parameters given in Table [T] 

analysis can be found in [llj . A detailed explanation of several advanced ideas is given in lOj. Our al- 
gorithm parallels one of the approaches from sensitivity analysis by using a Monte-Carlo type method. 
However, the key difference between our approach and known methods is that we take advantage of 
the vector field "directly" in computations to find a more efficient algorithm. In particular, we are not 
interested in a black-box scheme which works for all mathematical models, but to design a method 
which uses the structure of the model given by a continuous or discrete-time dynamical system. 



2 Limit Cycles in ODEs - The Spruce Budworm 

We use a classical example from mathematical ecology to illustrate our algorithm. The 2-dimensional 
ODE model proposed in [7] represents the interaction between the spruce budworm N and the leaves 
of a forest R. It is a simplification of a 3-dimensional system for the same ecological situation given in 
[B] . The paper [B] is partly based on earlier work in [5] , whereas re-expositions of the spruce budworm 
model can be found in many texbooks by now (see e.g. [13]). Our model reads: 
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R = piR[l ]-p7N 
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where R is the branch area and N is the budworm density. To keep the notation simple we deviate 
from the conventions in ecology and simply index the seven parameters pi, . . . ,pT. 



The descriptions, units and default values are summarized in Table [TJ We consider the model only 
for R,N>0 representing biologically sensible nonnegative quantities. 

Table 1: Parameters for the sample problem 



Parameter 


Description 


Units 


Default Value 


Pi 


intrinsic branch growth rate 


1/yr 


0.15 


P2 


intrinsic budworm growth rate 


1/yr 


1.6 


P3 


maximum branch density 


branches/acre 


24000 


P4 


maximum budworm density 


larvae/branch 


200 


P5 


maximum budworm predated 


larvae/ (acre*year) 


28000 


P6 


half of maximum density for predation 


larvae/branch 


1.5 


P7 


consumption rate of budworm 


(branches*acre)/(yr*larvae) 


0.0015 



We investigate the characteristics of a hyperbolic limit cycle 7 C M^, which represents long-term 
coexistence for branches (i.e. trees) and the budworm population. A stable limit cycle occurs at the 
given default parameters (see Figure [T|) . We observe from equation ([1]) and Table [T] that the variable 
R is slow compared to the variable N. This means that (H]) is a slow-fast system. Direct calculation 
of the nuUclines shows that the observed limit cycle is a relaxation oscillation. 

For equation ([1]) we have to explore a 7-dimensional parameter space. We denote the parameter 
space by P C M.l_ and an element in P by p = (pi, . . . ,P7). To indicate the dependence of 7 on p we 
shall also refer to the limit cycle as 7(p). 

Biologically relevant questions are to find the maximum values of the coordinates of 7(7)). These val- 
ues correspond to the maximum branch area Rmaxilip)) and maximum budworm density Nmaxilip)) 
along a cycle. The maximum budworm density is particularly important as it describes the maximum 
insect outbreak in the ecosystem. Our main focus is to find the directions in P which correspond 
to relevant variations of Rmaxilip)) or Nmax{l{p))- Notice that our aim is to determine parameter 
dependencies globally as well as locally, i.e. we are interested in the overall behaviour as well as sig- 
nificant variations around points. A further goal is to identify functional dependencies of Rmaxilip)) 
or Nmaxilip)) on subsets of parameters to distinguish between a linear and a nonlinear response. 

3 An Algorithm for Parameter Space Exploration 

We are going to describe the algorithm in the general setup of ODEs and illustrate all steps with the 
model given by equation ([T]). 



3.1 The Setup 

Consider an ordinary differential equation 

X = f{x,p) for a; G R", p e P C R'^' and / G C^{W+^) 



(2) 



where p is a vector of parameters p = {pi,p2, . . . ,pk)- Suppose that equation ([2]) has a quantitative 
property 7(p) such as a limit cycle (or e.g. a heteroclinic orbit). By "quantitative" we refer to the 
possibility of measuring an important characteristic of 7(p) like cycle amplitude (or e.g. length of a 



heteroclinic connection) . Furthermore assume that we have a compact set B C P in which the feature 
is persistent, i.e. we do not have any bifurcation points inside B. We shall refer to B as a hyperbolic 
domain for j{p). 

Define a function C : B —^ M. given by mapping a parameter vector p to C{p), where C{p) is the 
quantitative measure of j{p). Assume we have a method to compute C numerically. Notice that 
C might be very expensive to compute as a single function evaluation usually requires at least one 
complete numerical solution of equation ^ . 

In our sample problem we take B = {p-'^-'"'" * ± Si} C Rl_ where 6i are chosen based on the 
combinations and ranges of parameters we are interested in and to preserve the limit cycle j(p) found 
in the first quadrant for the default values given in Table [TJ If we compute "f{p) numerically we obtain 
a finite sequence of values {'y{p)R,k,7{p)N,k), where the indices "R" and "N" denote the coordinate 
directions in ([1]). If we are interested in the maximum insect outbreak we define C{p) by 

C{p) = max7(p)7v,fc (3) 

k 

Defining C{p) obviously depends on the given problem and the questions of interest in the mathemat- 
ical model. 

3.2 Discretization 

Assume without loss of generality that B is connected; if this is not the case we can compute for each 
of the components of B separately. Choose a finite grid G C B oi spacing Si in direction pi defined 
by intersecting a scaled lattice with B so that 

G = {(si • Z, S2 • Z, . . . , sfc • Z) e M*^} n B (4) 

Note that we have chosen the lattice to be uniformly spaced in each direction. This is not necessary 
for any of the following calculations. We should change to a nonuniform grid if there is a priori infor- 
mation about parameter regions available having large variations in C{p). In our model problem we 
have chosen a uniform spacing in each direction. 

Now consider the graph formed from G by connecting all points which differ in one coordinate pi 
by precisly Si. We also use G to denote this graph. Metrize G by 

d{p, q) = minimum number of edges connecting p to g (5) 

3.3 Relevance Measure 

The next step is to compare two given elements p,q ^ G and estimate when \C{p) — C{q)\ is "large". 
To do this we use the given vector field /. In the following we shall describe the most basic version of 
the procedure. 

Step 1: Select a set of points in phase space determined by 7(p). This is best illustrated by a 
concrete example. For our limit cycle 7 we can choose fci different parameter values. We compute 7 
for each parameter value and we select m points in phase space on the computed limit cycle. This 
gives a collection of points in phase space: 

Sx '■= {xi\ i ^ 1,2,.. . ,k2 ^ mki} (6) 



We choose fci and m such that the number of elements in Sx is small. In our model problem we 
have computed 7 for a set of values chosen uniformly at random from G with fci = 4. Then we 
selected to = 4 points on each computed limit cycle starting with a random point and uniform spacing 
along the cycle. This gives a sample representing the different shapes of possible cycles in phase space. 

Step 2: Choose a set points from parameter space 

Sp^{p\p',...,p''^}cG (7) 

uniformly at random. Fix a neighbourhood size in G, i.e. pick an integer Usize such that a neigh- 
bourhood of a point p is given by flp = {q Cz G\d{q,p) < Usize and q 7^ p}. For each value p^ d Sp 
choose a fixed number of neighbours ^4 at random from Qpj , call this set Spj ; usually we want to pick 
between k^ between 2 and 4, but any number less than the number of elements in flp is possible. In 
preparation for the next step, we mention that we should choose #(5*^) = fcs such that a summation 
of ki ■ k^ ■ ^2 elements is computationally feasible. 

Step 3: Estimate a measure of relevance. Again we consider the model problem for illustration. 
Let 
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(8) 



So r represents the average variation between two neighbours along the sampled points and hence 
can be taken as an estimator for the variation of the limit cycle from one point to a neighbouring 
point. Hence we have found a way to estimate when \C{p) — G{q)\ is "large". We call r the (average) 
relevance measure. We can also define 

M{p)^^Y^\\J{x,,p)\\, (9) 

which we refer to as the relevance function. We can easily compute the relevance function for different 
values of p as we only need a small number of evaluations of the vector field if Sx is small. We shall 
discuss more advanced methods of defining r and M{p) in Section ID 

3.4 Neighbour Comparison 

Let p,q ^ G be given in parameter space and suppose we have already computed G{p). Now we 
choose a neighborhood fip of p in G excluding p itself, e.g. Op = {q € G\d{p^ q) — 1}. Fix a tolerance 
tol S [0, 00), then the rationale for approximating the desired function G is: 

• If g G rjp and \M{p) — M{q)\ > tol ■ r ^ q shows a relevant variation so we compute G(g) 

• li q E Up and \M{p) — M{q)\ < tol ■ r ^ q does not show a relevant variation so set G{q) — C{p) 
Clearly we have two special cases: 

• tol = Q: In this case we always compute any point encountered by the algorithm exactly. 

• tol — l-\- maxp^ggG \M{p) — M{q)\: In this case we never compute any neighbours of points and 
just set their values to be the value of the center p. 

In our model problem we used the criterion for several values of tol between 0.9 and 1.1. The reason 
for comparing with neighbours will be explained in sections 13.61 and 13.71 



3.5 Iteration 

We describe next how to run the exploration algorithm to compute C{p). We have two choices to 
proceed: 

1. Fix a subset Go of G, for each point p G Gq compute C{p). Fix a neigbourhood size for each 
point such that neighbourhoods of points in Go do not intersect and perform neighbourhood 
comparisons and computations as described in section 13.41 We call this method deterministic 
exploration. 

2. Start with a random initial point p £ G. Compute C{p) and mark p as checked. Consider a 
neighbourhood Op and perform neighbourhood comparisons and computations as described in 
section [3^ for all neighbours which have not been marked yet. Mark all elements in Op. Now 
choose a new random point out of all unmarked points. We call this method (uniform) random 
exploration. 

For both methods we can choose a fixed number of points imax, which should be explored. Clearly 
imax should be defined so that imax < #(G). Having reached imax we have obtained the function G 
on a subset G* of the grid G. Notice that if the grid size shrinks to zero and imax approaches #(G) 
the algorithm converges to a full computation of G. 

The most successful algorithms for higher-dimensional numerical problems often use an element 
of randomization. A classical example is numerical integration. In this area Monte-Carlo techniques 
are predominant for any problems of dimension > 4 (see e.g. 2 for a survey). Hence we expect the 
deterministic exploration to work only for small number of dimensions and recommend the random 
exploration. 

3.6 Interpolation 

Given the values of G on G* we obtain all values in the hyperbolic domain B by interpolation. This 
constructs G on the full domain B. For our model problem we used linear interpolation. For details 
on basic numerical methods for interpolation consider ^T . Most algorithms currently in use for multi- 
dimensional interpolation are based on a Delauney triangulation of given data. Many use the 'quick 
huU'-algorithm (QHULL, see [I],[H]); we used the version of QHULL implemented in MatLab for our 
sample problem. 



One reason why we have used neighbour comparison as described in section 13.41 is obvious from 
the one-dimensional toy-example illustrated in Figure [21 Independent of all previous considerations, 
just suppose the 1-dimensional interpolation problems on [0, 1] are: 

{(0, 0.3), (0.5, 0.25), (1, 0.3)} Problem 1 
{(0, 0.3), (0.5, 0.25), (0.58, 0.04), (1, 0.3)} Problem 2 

One extra point gives the indication that there is a relevant decrease of the function value to the 
right of 0.5 whereas no neighbour comparison would entirely miss this decrease. This means that the 
local comparion helps us to detect globally important directions of variation of j{p), which could be 
missed just using point-wise Monte-Carlo exploration without neighbour comparison. It also indicates 
directions which are of interest near a point. We refer to this feature of our method as local/global 
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Figure 2: Interpolation for 3 points linear, 4 points linear and 4 points cubic spline 



approximation. 

Furthermore we see that other interpolation methods, such as the cubic spline in Figure [21 work 
as well, but might introduce a functional dependence due to the interpolation method. Therefore the 
conservative approach is to use linear interpolation and investigate interesting regions in parameter 
space more closely once a significant decrease or increase in the function C is found. 

3.7 Results and Visualization 

To determine the influence of parameters on the structurally stable feature 7(p) we have to organize 
the output of the previously described steps, i.e. we can either work with a subset of values G* or its 
interpolation to the hyperbolic domain B. Some problems arising are: 

• Arc there "enough" points to perform interpolation on B or on some given subset of Bl 

• How to visualize the output of high-dimensional data? 

We shall not discuss the two questions in detail as they form subjects on their own and just focus 
on a few special issues. 

The method of neighbour comparison balances the effect of dimension with respect to the number 
of points necessary for interpolation. Consider our model problem and suppose we fix five parameters 
on G* and want to interpolate C{p) over the remaining two parameters, i.e. over a 2-dimensional 
"slice" of B. Using just Monte-Carlo simulation without neighbour comparison will inevitably lead 
to many slices, which do not have any points on them. With neighbour comparison, by possibly in- 
creasing the neighbourhood size, we can assure that we have enough points to perform interpolation. 
Note that higher-dimensional interpolation techniques can also fail if there are too few points. For 
example, a Delauney triangulation might not be supported by many algorithms for sparse G* (see [T], 
[8] and www.qhull.org). 

For a fc-dimcnsional parameter space with G* C K*^ we can fix A; — 2 parameters and plot 3- 
dimensional output for C{p). This can be achieved by viewing G{p) as a surface over the two variable 
parameters. We found this method particularly useful coupled to animations, which allow for moving 
a third parameter and hence produce moving surfaces. During the algorithm we can easily keep track 



of points, which have been computed and those, which have been set to a local constant by neighbour 
comparison. Therefore we are able to identify regions and directions containing relevant variations. 

It should be clear that it is also easy to get estimates of maximum and minimum values over G* 
and hence over B. For example, in the budworm problem posed in section [5] we can now determine 
for any given region of parameter space an estimate of the maximum insect outbreak. 

3.8 Summary of the Algorithm 

The main steps can be summarized as follows: 

1. Identify the quantitative feature 7(p) and its hyperbolic parameter domain B. 

2. Define a quantitative measure of 7(73) encoded in a function C : S ^ M. 

3. Discretize the space B and define a metric on the discretization G. 

4. Sample points from phase space "representing" the feature 7(p). 

5. Choose a measure (e.g. the norm) to estimate the variation of 7(p) from the vector field /. 

6. Compute the relevance measure r and setup the relevance function. 

7. Choose discrete or random exploration and a size of neighbourhoods for computed points. 

8. Define a neighbour comparison based on r and the relevance function. 

9. Iteratively compute C on a certain number of points using exploration. 

10. Interpolate (linearly) between the values obtained for G . 

11. (Visualize and/or evaluate the interpolated function G to explore parameter dependence.) 

The overview shows that the algorithm works for other dynamical systems like maps, difference 
equations or partial differential equations. 



4 Advanced Local Exploration 



Instead of using the norm of the vector field, we can also use derivatives in the parameter directions, 
which leads to the relevance measure: 
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(10) 



where we recall that k is the dimension of the parameter space. The corresponding relevance function 
is given by: 



1 



^'^ipy-= TnT.1: 



=1 1=1 



df, . 

Pi 



(11) 



Remark: Using the partial derivatives is a technique employed in local sensitivity analysis and has 
the potential to produce a better measure of relevance. We can use finite-difference approximations 
to the derivatives if there is no explicit formula of f{x,p) avaiable. Mq{p) has the drawback that the 
computational cost increases. 



a u.t 




' 






c 










a> 










E „„^ 




-* 






■■= 0.35 










CO 




















c 










o 










it "-3 




* 


* 


■ 












o 






























en 










d 0.25 








* 














# 


















Q- 0.2 


-*r 






. 


E 










o 










CO 










O 











o 

" 4000 






* 


0) 

E 








c 3000 
o 






■ 


15 
3 
§- 2000 




* 
* 




o 
o 








^ 1000 

3 


# 




■ 



0.5 



1 1.5 

grid-points 



2 2.5 

xio" 



0.5 



1 1.5 

grid-points 



2 2.5 

xio" 



Figure 3: Comparison of computing time for expforation algorithm and full caleulation 

5 Numerical Results 

Note that since equation ll]) is a fast-slow system, it is stiff and explicit numerical solvers are going to 
fail. We used a modified Rosenbrock method of order two with explicitly supplied Jacobian matrix; 
see [9] for historical reference, [12] for the implementation and [4j for a modern exposition of the 
Rosenbrock method. We have simulated the system for the following range of parameters: 



Table 2: Ranges of parameters for the sample problem - spruce budworm model 



Parameter 


Default Value 


Range 


Spacing Si 


Pi 


0.15 


[0.149,0.151] 


0.0001 


P2 


1.6 


[1.5,1.7] 


0.1 


P3 


24000 


[22000,26000] 


250 


Pi 


200 


[190,210] 


2 


P5 


28000 


[24000,32000] 


500 


Pa 


1.5 


[1,2] 


0.1 


P7 


0.0015 


[0.001,0.002 


0.0001 



The main numerical question to pose is, what computational cost do we save by using our algo- 
rithm? This is hard to answer precisely as our goal is to provide a method exploring parameter space, 
which allows applied researchers to detect the main quantitative dependencies in their models. To 
perform a practical test we have simulated 10 subsets of the full parameter space grid defined in Table 
[21 The results of this simulation are shown in Figure [3] We have chosen to use random exploration of 
five percent of the total number of grid points; together with neighbourhood exploration this yields a 
sufficiently dense set of values to investigate the function C . 

From Figure[3|we observe that the computing time for our algorithm lies between 20 and 40 seconds 
whereas a full grid computation can take up to 4800 seconds. Now we have to answer the question 
if we can really extract all the quntitative information we are interested in from our algorithm in 




Figure 4: Maximum budworm outbreak as a function of p^ and pg for ps = 24000, 28000, 32000 - 
algorithm 

comparison to the much slower full grid computation. 

Suppose we are interested in the relations between p3, p5 and pQ. We simulate the system with 
P3,P5 and pq in the ranges given in Table [2] and use default values for all other parameters as given 
in Table [TJ We plot the maximum budworm outbreak N^ax as a function of p^ and pe- Plots for 
P5 = 24000,28000,32000 are shown in Figure [H The red dots indicate grid points which obtained 
their values from neighbours. The black dots indicate points which have been fully computed. The 
computation has been performed with a relevance measure r given in equation ([H]) and tolerance 
tol — 1.1 (see section [ 



Just given Figure [5] we can draw several conclusions: Obviously the outbreak is decreasing with 
increasing value oi pe. It is not significantly affected by variations of p3 compared to p^ and pe. We 
can see that increasing p^ decreases the maximum outbreak by comparing among all three figures. 
The plots show a weakly nonlinear relationship between pa and pe. We can also conclude that the 
quantitative feature varies "uniformly" in all three parameter directions, i.e. does not exhibit several 
nonlinear relationships in our region of interest. 

Now let us compare these findings with the computation on the full grid, which is given in Figure 
[5l We see that all our previous conclusions have been verified by the full computation. The only obser- 
vation we should make is that the first plot in Figure |3] and Figure [5] do not seem to match well. This 
results from the fact that our algorithm has produced few values for small pe and small values of ps 
so that this region was not interpolated. We could simply remedy this by extending the current inter- 
polation to the boundary but we decided on the most conservative approach given the computed data. 

Comparing the computational cost we have computed the two figures on several standard desktop 
workstations and obtained that it took approximately 10-15 seconds to compute with our algorithm 
to obtain NmaxiPsTPbTPe) for the hyperbolic domain B given by the ranges in Table [51 This means we 
have obtained a local/global approximation to N^ax on 32 • 16 • 10 = 5120 points. The computation 
of a single slice of pa and pg for a fixed value of ps took approximately 45-50 seconds on the same 
machines, yielding a total time for the same amount of points given by > 45 • 16 = 720 seconds. This 
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Figure 5: Maximum budworm outbreak as a function of p^ and pe for p^ — 24000, 28000, 32000 - full 
computation 



means that in this practical example we have used only a fraction of the full computation 



15 
720 



1 

48 



(12) 



This 50-fold increase in computation time supplied us with the same conclusions we would have ob- 
tained by "brute force" . Notice that this effect increases with dimension. 

The computation of a single slice of p^ and pg using our algorithm is about 1 second on most 
current workstations. This means that we can "move" with 3-dimensional plots through parameter 
space in real time and adjust the parameters along directions of interest according to the given output. 

We have also performed calculations for the full 7-dimensional parameter space and obtained 
local/global approximations of Nmax and Rmax ■ Although this computation takes around 2 hours for 
most values of tol Ki 1, it is still feasible. We remark that the algorithm can also be 'scaled' according 
to the computational tools available and the desired goals: 

• Accurate computations, high computing power, long-waiting time: We decrease tol and select a 
relatively large fraction of sample points Sx ■ 

• Coarse global estimates, low computing power, direct feedback: We increase tol and reduce the 
size of Sx ■ 

This means that powerful computing environments can be very useful, but are not necessary to extract 
reliable information from a model. 

6 Error Analysis 

The precise error analysis is not the focus of this paper. Nevertheless we shall motivate in this section 
why we expect an error that behaves like 1/Vn where n is the number of computed gridpoints. 
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Assume that the hyperboHc domain B is a compact hyperrectangle. First consider the discretizations 



Gj C B C M*^ such that the foUowing condition holds: 



min \\p — q\\2 -^ for all g G _B as j ^ oo (13) 

peGj 

where ||.||2 denotes the Euclidean norm in K.'^. We call the sequence of grids Gj an approximation 
sequence for B if condition (fTS]) holds. li f : B ^ M. then we define the usual L^ norm as: 



II/IIli = / \fip)\dp (14) 

Jb 

Since i? is a compact hyperbolic domain we conclude that C : B — > M is bounded (for the definition 
of the function C, see l3.ip . Let Ca : B -^ R denote the local/global approximation to G obtained by 
our algorithm for the grid Gj. Recall that we have computed n{j) points (depending on imax and the 
neighboourhood size) on the grid and used linear interpolation to construct Ga on B; denote these 
points by {a^fcl^^i- Obviously Ga is a random variable depending on the relevance measure r and the 
grid Gj. Hence we should formally write Ga — GA{r.j), but for simplicity of notation we sometimes 
just write Ga. Since G is smooth we can modify Ga on the complement of Gj to be piecewise linear 
and satisfy Ga{p) > G{p) for all p & B. 

Let P denote the probability measure on B induced by a uniformly distributed random variable 
on B. Since B is compact we can assume without loss of generality that J„dx = l so that the density 
of P is identically 1. Now define 

-, "(i) 

^■^' fe=l 

We also set I — J^ G{p)dp. This allows us to consider a standard result from Monte-Carlo integration. 

Proposition 6.1. If tol = then GAir.j) = Ga(j), *-e- Ga is indepedent on the random variable r. 
Suppose n{j) -^ oo as j ^ oo then 



^AU) - I 



<^/VnU) 



as j -^ CO (16) 



where Z ^ A^(0, 1), i.e. we have convergence in distribution to a normal distribution with mean and 
variance 1. 

Remark: Note that j —^ oo implies n{j) -^ oo even if we only increase the number imax linearly. 

Proof. (Sketch, see [5]) If tol = 0, no neighbour comparison is performed and all elements are com- 
puted, so that Ga is independent of r giving the first result. 

Notice that Var{lA(j)) is bounded since G is bounded. Then the central limit theorem implies the 
desired result. D 

For two sequences of random variables {Xn}, {F„} we write X„ = Op{Yn) if for any e > there 
exist Me > 0, A^e > such that 



P 



Xn 
Yn 



> M,] <e for all n>N^ (17) 
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so that Xn/Yn is bounded in probability. Therefore we obtain from Proposition 16.11 that : 
Since Op{.) ~ 0{.) and Ca{p) > C{p) for all p e B we can write 

|..o,-^IH|C.-qu.^o,(-^).o(-^) (19) 

where the last big-oh notation refers to the deterministic asymptotic behaviour on B. We can state 
the previous discussion as a practical result: 



Proposition 6.2. If tol = then \\Ca — C\\li ~ O I , j for n{j) (respectively j) sufficiently 

large. 

The convergence of the error by an inverse square-root law is obviously very slow. Nevertheless, 
it is precisely the law we would expect from the use of Monte-Carlo methods. It reflects that we are 
dealing with a problem without restriction on dimension. Proposition 16.21 justifies our initial claim 
that we expect an error that behaves like 0(-^). 

7 Conclusions 

We presented a new algorithmic paradigm to explore the quantitative variation of structurally stable 
features in dynamical systems depending on parameters. We illustrated our algorithm in the case of 
ODEs and a 2-dimensional model problem and discussed the main steps of discretization, sampling, 
relevance measures, neighbour comparison, random exploration and interpolation. Numerical results 
demonstrate the practical value of our method and its efficiency compared to direct computations. 
Furthermore, error estimates for the random exploration strategy were proven and confirmed the ex- 
pected Monte-Carlo type behaviour with error of order 0{l/y^) independent of the dimension of the 
phase space or parameter space. 

Despite the detailed error analysis, we stress that our primary goal was not to generate accurate 
data but to provide practioners with a computationally efficient method. The efficiency derives from 
a strategy which has been very successful in the analytical description of dynamical systems - we do 
not compute an explicit solution of an equation but characterize its properties directly from the equa- 
tions and their structure. We have shown that this general idea can be very useful for the numerical 
investigation of dynamical systems too. 
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